CalcSparesRisk<-function(BNcK,PartProbs,RecertSched,PartCost,ThroughPut,Capacity,RecertCost,StorageCost,alpha,scale) {
  
  
  # This function is designed to calculate the disparity between purchased spares and demand for spares 
  # for a given system based on part reliability, recertification schedule, part indenture in the BOM, and purchase plan
  # Function starts at year 1 and adds up the current spares for each part, subtracts parts
  # the at predicted to be used, adds those parts to be canniballized from replaced parts (Harvesting is not currently implemented in the code.
  # It then outputs the total number of spares to be purchased each year to maintain an inventory suitable to
  # meet the no stock-out confidence facotr desired (Demand).  This Demand is compared to the actual purchae plan to create a cost impact of
  # under or over buying parts
  #   
  # The inventory level is determined thorugh use of a CASA program equation that assumes an
  # exponential failure rate (as this program does in the RPP) and determines the inventory amount
  # to meet the spares demand of this failure rate with a given confidence.
  #      S : x = 0->s [Σ((e^-d)(d^x)/(x!)) >= c]
  #   Where d = failure rate, c = confidence, s = inventory stock
  #
  # It is assumed that the inventory begins with two of every part (program plan to purchase
  # two entire missiles as spares) 
  #   
  # Number of each part replaced is calculated with the probabiltiy of replaced from BNcK (branched
  # probability of replacement matrix) and the number of missiles to be recetified.  this iteration
  # assumes no distribution of the part replacement probabilities.
  # 
  # Cannibolized children are calculated first by identifying the children of a replaced part,
  # then finding if any of the children also fail (using the PartProbs value), and adding the 
  # surviving children to the inventory.  Cannibolized parts have a shelf life of less than a year
  # and are removed from next year's inventory.
  # 
  # It is assumed that all spare parts have a shelf life of 4 years, the program will remove the 
  # initial spares from three years prior (the final spares from 4 years prior)
  # 
  # Inputs
  #   PartTree,(not implemented yet )
  #   PartProbs,
  #   BNcK,
  #   ParentIdx,(not implemented yet )
  #   ChildIdx,(not implemented yet )
  #   Conf
  #   Recert
  # 
  #
  # Output
  #   CostImpact -  A matrix that contains to the cost of overbuying or underbuying spares.  Overbuying cost is the total cost
  #                 of unused parts that are purchased according to the current plan.  Underbuying cost is the cost associated with
  #                 wasted time awaiting a missile (addition cost for ship in port, cost of missile in inventory, wasted parts on missile )
  #          
  #
  # Written by: Nicholas J. Bouchat, Willcor Inc., December 2012
  #___________________________________________________________________________________
  
  #Initialize Variables
  RecertSched = read.csv("RecertificationTotals.csv",header=FALSE, sep=",", quote="", dec=".",as.is=TRUE) #read in csv build from recertification schedule
  Recert<-as.integer(RecertSched) #transform Recert Schedule data to a vector of integers
  #may want to change this to the delivery schedule and mod delivery schedule to recert schedule within the code
  
  PartCostI = read.csv("PartCost.csv",header=FALSE, sep=",", dec=".",,,)  #read in part cost matrix from CSV file
  PartCost <-as.numeric(PartCostI) #transform Cost data to a numeric vector
  PartCost<-PartCost[-1]   #truncate matrix to remove missile top level
  PartCost<-PartCost[-1]   #truncate matrix to remove first stage level
  
  PartProbsI = read.csv("PartProbs.csv",header=FALSE, sep=",", dec=".",,,)#read in part reliabilities from csv file
  Reliability <-as.numeric(PartProbsI)  #transform reliability data to a numeric vector
  PartProbs = 1 - Reliability  #define the failure rate as 1 - Reliability
  PartProbs<-PartProbs[-1]  #truncate matrix to remove missile top level
  PartProbs<-PartProbs[-1]  #truncate matrix to remove first stage level

  load("SM3_Demand.RData")
  
  N = length(PartProbs)  #identify the total number of parts for use in the loop
  OpYears = length(Recert) #identify the number of years to run by using the recert schedule
  RecertCost = 69000 #Cost per missile for recertification at Huntsville (provided by program)
  Throughput = 1  #throughput of the recertification facility (amount of time it takes to recertify a missile in months)
  Capacity = 5  # the number of missiles that can be recertified at once
  StorageCost = 1000000 #cost per missile per month for missiles to be stored at the facility
  ShelfLife=4   #the number of years parts can sit on the shelf before they an unusable in a feilded system
      
  BNcKProb <-BNcK$Prc #initialize the Demand probability to be equal to the part reliability
  BNcKProb<-BNcKProb[-1]  #truncate matrix to remove missile top level
  BNcKProb<-BNcKProb[-1]  #truncate matrix to remove first stage level
  
  Demand = roundUp(outer(BNcKProb,Recert),1)  #simple calculation of number replaced rounded to next integer.  Can change this to distributed method
                                  #N by t vector has all predicted replacements for a given part along a single row
  
  Purchase=CalcPurchase(alpha,scale,BNcKProb,Recert,Demand)
  
#   Purchase <- matrix(data=NA,nrow=N,ncol=OpYears)
#   Purchase = outer(PartProbs,Recert)  #this is the planned purchase rate matrix (buy each part as its probability of failure dictattes) This will be defined by the user

  RecertRate = Capacity/Throughput*12  #calculatte the yearly throughput for recertification.  This can be refined to a monthly rate for better fidelity
  
  ExpiredSpares = CalcExpiredSpares(OpYears,N,Purchase,Demand)
      
  Stock <- CalcStock(OpYears, N, Recert, Purchase, Demand,ShelfLife)
  
  Conf <- CalcConf(OpYears, N, Stock, Demand)
      
  Diff = Stock - Demand    #Diff is calculated as the difference between the supply of spare parts and the demnad for spare parts.  The need for parts tranlates to missiles that cannot be returned to the fleet
     
  MissileStock <- CalcMissileStock(OpYears, N, Diff, RecertRate, Recert)
      
  Cost <- CalcCost(OpYears, N, PartCost, Purchase, StorageCost, MissileStock, Diff, RecertRate, RecertCost)

#   Risk = list(Conf=Conf,Cost=Cost)
  Risk = c(Cost,Conf)
  
  return(Risk)
}#end spares risk function

#=================================================================================================================================
CalcPurchase<-function(alpha,Scale,BNcKProb,Recert,Demand){
  
  Parts = length(BNcKProb)
  Year = length(Recert)
  Avg = median(BNcKProb)
  
  InvProbs = max(BNcKProb)+min(BNcKProb)-BNcKProb
  
  PurchaseRate = alpha*BNcKProb+(1-alpha)*InvProbs
  
  Purchase = roundUp(outer(PurchaseRate,Recert) + Scale,1)
  BOMPart <- seq(1,Parts,by=1)
#   for(year in 1:Year){
#     NumParts = Purchase[,year] 
#     Demands = Demand[,year]
#     plot(BOMPart,NumParts,col= "blue", ylim=c(0,10))
#     lines(BOMPart,NumParts,col="blue")
#     lines(BOMPart,Demands,col="red",pch=22,type="o")
#     title(main = year)
#     legend(75,10,c("Purchase","Demand"),cex=.8,col=c("blue","red"),pch=21:22,lty=1:2)
#   }
  
  return(Purchase)
}
#========================================================================================================
#This function is currently hardcoded for a shelf life of 4
CalcExpiredSpares<-function(OpYears,N,Purchase,Demand){
  
  ExpiredSpares<-matrix(nrow=N,ncol=OpYears)
  LeftoverStock<-matrix(nrow=N,ncol=OpYears)
  Expire<-matrix(nrow=N,ncol=OpYears)
  
  for(Part in 1:N){
    for(t in 1:3){
      ExpiredSpares[Part,t]=0
    }
    for(t in 4:OpYears){
     
        if (t == OpYears){
          LeftoverStock[Part,t] = sum(Purchase[Part,])-sum(Demand[Part,])
          FinalExpire[Part,t] = Purchase[Part,t-3]-Demand[Part,t-3]-Demand[Part,t-2]-Demand[Part,t-1]-Demand[Part,t]
          ExpiredSpares[Part,t] = FinalExpire[Part,t]*(FinalExpire>0)+(LeftoverStock[Part,]*(LeftoverStock[Part,]>0))
        }
        else{
          
          Expire[Part,t] = Purchase[Part,t-3]-Demand[Part,t-3]-Demand[Part,t-2]-Demand[Part,t-1]-Demand[Part,t]

        }
     }
 }
  ExpiredSpares<-Expire*(Expire>0)
  return(ExpiredSpares)
}

#====================================================================================

#Stock function to determine how many missiles are in Stock for a given year.  Based on the Demand, purchase schedule, and part shelf life.
CalcStock<-function(OpYears, N, Recert, Purchase, Demand,ShelfLife,ExpiredSpares){
  
  Stock <-matrix(nrow=N,ncol=OpYears)  #initialize the vecotr to contain the number of parts in inventory (Stock)
  
  Stock[,1] = 0     #set initial stock value to 2
  
  
  
  for (t in 1:(OpYears-1)) {
    for(Part in 1:N){
          
#       cat("StockSize",dim(Stock),"\n")
#       cat("Stock",Stock[Part,t],"\n")
#       cat("Demand",roundUp(Demand[Part,t],1),"\n")
#       cat("Purchase",roundUp(Purchase[Part,t],1),"\n")
#       cat("t",t,"\n")
#       cat("Part",Part,"\n")
            
      if (t < ShelfLife){   #run the spares purchase calculation for the first x years 
        #during which time no spares should reach shelf life
#         cat("Within Shelf LIfe \n")
        Stock[Part,t+1] = Stock[Part,t]-Demand[Part,t]+Purchase[Part,t]
      }
      else {    #run the spares purchase calculation for the number of years in recert vector
#         cat("Outside Shelf LIfe \n")
                
        Stock[Part,t+1] = Stock[Part,t]-Demand[Part,t]+Purchase[Part,t]-(ExpiredSpares[Part,t]*(ExpiredSpares[Part,t]>0))
        
      }
          }
    BOMPart<-seq(1,N,by=1)
    NumParts = Purchase[,t] 
    Demands = Demand[,t]
    InStock = Stock[,t]
    plot(BOMPart,NumParts,col= "blue", ylim=c(0,10))
    lines(BOMPart,NumParts,col="blue")
    lines(BOMPart,Demands,col="red",pch=22,type="o")
    lines(BOMPart,InStock,col="green",pch=23,type="o")
    title(main = t)
    legend(75,10,c("Purchase","Demand"),cex=.85,col=c("blue","red"),pch=21:22,lty=1:2)
    #above is code to print the demand and purchase curves when alpha and scale are not implemented.
#     NegStock=Stock*(Stock<0)
#     StockNegs<-NegStock[,t]
#     BOMPart<-seq(1,N,by=1)
#     plot(BOMPart,StockNegs,ylim=c(-12,0))
#     lines(BOMPart,StockNegs)
#     title(main=t)
  }
  
#   cat("Demand",roundUp(Demand,1),"\n")
#   cat("stock",Stock,"\n")
  return(Stock)
}
#end Stock

#====================================================================================

#confidence function calculates the "no stock out confidence" based on the calculated demand and current spares inventory
CalcConf<-function(OpYears, N, Stock, Demand) {
  
  PartConf = matrix(nrow=N,ncol=OpYears) #initialize part confidence matrix to store part stock confidence for each part for each year
  Conf = numeric(OpYears)
  Stock4Conf<-Stock    #initialize stock matrix to be used for confidence calculation
#   DPt=roundUp(Demand,1)
  
  for (t in 1:OpYears) {
    for(Part in 1:N){
      
      Demand[Part,t]->DPt
#       if (is.na(DPt)){
#         print(Demand)
#         stop ()
#       }
      Temp= 0 #initialize a temporary vector to house values for confidence equation before performing summation
    
      if (Stock[Part,t]<0){#add line here to normalize negative stock values to zero.  Purchase cost is included in the underbuy cost (line 236)
        Stock4Conf[Part,t]=0
        }
      for (x in 0:Stock4Conf[Part,t]){
#         cat("t",t,"\n")
#         cat("Part",Part,"\n")
#         cat("x",x,"\n")
        
        Temp = Temp + (exp(-DPt)*DPt^x)/factorial(x)
        
      }
      
      if (Temp == 0){ #error is being caused by negative stock values being fed into the factorial above
        Temp=1
      }
      
      PartConf[Part,t] = Temp
    }
    Conf[t] = min(PartConf[,t])
    
  }
  
  TotalConf = min(Conf)
  return(TotalConf)
}
#end Conf

#========================================================================================

roundUp <- function(x,to) {    #function used to round up valudes in the Demand and and purchase matricies to integers

  y = to*(x%/%to + as.logical(x%%to))
  
  return(y)
}

#========================================================================================

#Missile Stock takes into account insufficient recert rate and insufficient spare parts to calculate the number of missiles that will be unable to be returned to the fleet
CalcMissileStock<-function(OpYears, N, Diff, RecertRate, Recert) {
  
  MissileStock = numeric(OpYears) #initialize vector to contain the number of missiles in stock for a given year
  UnfixedMissiles = numeric(OpYears)  #initialize vector to contain the number of unfixed missiles due to insufficient spares for a given year
  
  for (t in 1:OpYears) {
    for(Part in 1:N){
      
      UnfixedMissiles[t] = max(abs(Diff)*(Diff<0)) #number of unfixed missiles due to insufficient spares calculated as the max number of parts needed for any part
                               
      if (RecertRate>Recert[t]) {  #find the number of missiles remaining in stock if Throughput and Capacity are sufficient but spares stock is not
       MissileStock[t]<-UnfixedMissiles[t]
      }
      #end insufficient stock stock w/ sufficient throughput & capacity
      else { #throughput and capcity insufficient
       MissileStock[t] = Recert[t]+UnfixedMissiles[t]-RecertRate
      }
    }
  }
  #end insufficient throughput and capacity
  return(MissileStock)
                           
}#end MissileStock

#========================================================================================

#Cost is a function that determines the cost associated with the acquisition of spares in the purchase plan and the recertification costs for the current plan 
CalcCost<-function(OpYears, N, PartCost, Purchase, StorageCost, MissileStock, Diff,RecertRate, RecertCost){
  
  AnnualOverbuyCost = matrix(nrow=N,ncol=OpYears)
  AnnualUnderbuyCost = matrix(nrow=N,ncol=OpYears)
  AnnualPartCost = matrix(nrow=N,ncol=OpYears)
  AnnualStorageCost = numeric(OpYears)
  AnnualFixedRecertRateCost = numeric(OpYears)
  AnnualCost = numeric(OpYears)
  Cost = 0 
  
  for (t in 1:OpYears) {
    for(Part in 1:N){
#       cat("t",t,"\n")
#       cat("Part",Part,"\n")
      AnnualOverbuyCost[Part,t] = ExpiredSpares[Part,t]*PartCost[Part] #Cost impact is wasted money used to buy parts that are not needed
      if(t>2){
      AnnualUnderbuyCost[Part,t] = abs(Diff[Part,t-2])*(Diff[Part,t-2]<0)*PartCost[Part] #cost to buy parts that were not purchased in sufficient quantity
      #-------------------------------------------------------^add an expedite factor to the cost if needed
      }
      else{
        AnnualUnderbuycost[Part,t]=0
      }
      AnnualPartCost[Part,t] = Purchase[Part,t]*PartCost[Part]  #annual cost of initial spares purchase
      
      AnnualStorageCost[t] = MissileStock[t]*StorageCost #annual cost to store missiles due to insufficient throughput or spares acquisition
      
      AnnualFixedRecertRateCost[t] = RecertRate*RecertCost*.2 #calcuate the cost of a given throughput and capacity as 20% of the recertification cost per missile
      
      AnnualCost[t] = sum(AnnualOverbuyCost[,t])+sum(AnnualUnderbuyCost[,t])+sum(AnnualPartCost[,t])+AnnualStorageCost[t]+AnnualFixedRecertRateCost[t] #Calculate the annual cost due to under/overbuying spares
      
      Cost = sum(AnnualCost)
    }
    
  }
  return(Cost)
}#end Cost function